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Legendary mathematician Paul Erdos liked to talk about the perfect proofs of mathematical theorems maintained 
T— I by God. In [1 1, six different proofs are given to the infinity of prime numbers. Unlike mathematical proofs which are 
O mentally reproducible objects, signal processing researchers can have an intimate discussion with God by creating 
^ computational algorithms which are experimentally reproducible objects. This article presents a class of projection- 
^ ! based solution algorithms to the problem considered in the pioneering work on compressed sensing [21 - perfect 
2 reconstruction of a phantom image from 22 radial lines in the frequency domain. Under the framework of projection- 
based image reconstruction, we will show experimentally that several old and new tools of nonlinear filtering all 
■ lead to perfect reconstruction of the phantom image, which suggests that the result achieved by /i -optimization is 
Q ■ less like magic. 



I. Background 
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In in, the authors reported a "puzzling numerical experiment" which obtained the perfect reconstruction of a 

m ■ 

O phantom image at the sampling rate of 50 times smaller than the Nyquist ratefj. This experiment has motivated 
^ the authors of [2| to obtain a nonlinear generalization of Shannon's sampling theorem [3| and others to develop 
computationally efficient algorithms for ^i-optimization. The rapidly increasing interest in "compressed sensing" 
(CS) has a significant impact on the community of signal processing - e.g., various special sessions and special 
issues have been devoted to this topic. At SPIE Conf. on Visual Comm. and Image Proc.'2010, Prof. Changwen 
Chen at SUNY Buffalo asked an insightful question to the panelists (including the author of this column paper): 
what new insight do you think CS - the hot computational tool - to signal processing research? This article is 
based on the follow-up thoughts about Prof. Chen's question; while the author has chosen to indirectly answer his 
question from an algorithmic perspective - i.e., by demonstrating how the problem of phantom image reconstruction 
can be solved by a number of standard tools including Perona-Malik (PM) diffusion f4], nonlinear diffusion lH, 
translation-invariant (TI) thresholding 16] and Shape-adaptive DCT (SA-DCT) |7]. In particular, it will be shown 
that the classical framework of alternating projection manifests as much magic as /i -optimization (if not more). 

'strictly speaking, the term of Nyquist rate is only defined for analog signals (not for discrete signals) in the literature of signal processing. 
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II. Nonlinear Image Filters as Nonexpansive Maps 

The limitation of linear filtering on image signals has been recognized as early as 1970s [8|. In late 1980s 
and early 1990s, two lines of attacks on nonlinear filtering became influential - nonlinear diffusion and nonlinear 
thresholding. The idea of nonlinear diffusion originated from the scale-space analysis of image signals developed 
by vision community ||9l - the key insight behind Perona and Malik's scheme f4\ was to introduce a set of nonlinear 
diffusion coefficients with edge-stopping capability. Even though the stability of Perona-Malik diffusion has mostly 
been shown experimentally, it has sparkled a whole new school of thoughts on developing nonlinear tools for image 
analysis and processing. Two years later, the model of total-variation (TV) diffusion was proposed by Rudin and 
Osher fTOl - the key insight in their work is to replace I2 by li. Another two years later, wavelet shrinkage or 
thresholding was proposed by Donoho 111 II - the key insight there is to recognize the importance of singularities 
or tails (of heavy tail distributions). The connection between wavelet-shrinkage and TV-diffusion was established 
much later in 2004 |[T2l . Similar ideas or extensions related to li have been rediscovered several times including 
least absolute shrinkage and selection operator (LASSO) (13], basis pursuit (BP) [14] and the latest compressed 
sensing (CS) d. 

Both diffusion-based and thresholding-based nonlinear filters can be viewed as nonexpansive maps ifTSl - i.e., 
\\Pf\ \ < I I/I I- Rigorous proof for the nonexpansiveness of thresholding operators is relatively easy to obtain (e.g., 
|[T6l ): while similar success has not been achieved for nonlinear diffusion operators (e.g., the convergence of Perona- 
Malik diffusion has shown to be notoriously difficult to establish analytically, which gives rise to so-called Perona- 
Malik paradox IITtI ). Similar sentimental comments can also be made about projection-based texture-synthesis lITSl 
where authors have found experimentally their algorithms converge for all test images despite the lack of convexity 
of the defined constraint sets. In the mathematical literature, Ekeland seems to be the first one recognizing the 
importance of nonconvex minimization; at the end of [19], he advocated "to seek some kind of saddle point instead 
of a minimum". The implication of this proposal into signal processing research seems to be: instead of articulating 
the objective function for minimization, one could design a constraint set (not necessarily convex) or its associated 
projection operator (even based on heuristics). Such projection-first point of view (in contrast to variational or 
energy-first) arguably better fit the taste of electrical engineers and computational scientists. 

III. Image Reconstruction via Alternating Projections 

The power of image reconstruction via alternating projections was discovered as early as 1978 by Youla [20|. 
The idea is extremely simple - as long as we can come up with more than one constraint set for a target, alternating 
projections onto constraint sets offers a plausible solution to approximate the unknown target. Under the context 



of image reconstruction, we often face two kinds of constraint sets: one defined by the observation data (e.g., the 
Fourier coefficients along 22 radial lines) and the other associated with image prior (a.k.a. regularization functional). 
Both diffusion-based and thresholding-based nonlinear filters can be interpreted as projection operators onto the 
prior constraint set; the subtle difference between them is often more on the transient behavior (e.g., the speed of 
convergence and the route toward it) than their asymptotic one (since both approaches reflect the a priori knowledge 
about the signals of our interest). A common trick for improving the numerical stability of alternative projection- 
based algorithms P is called deterministic annealing - i.e., one can gradually decrease the diffusion coefficient or 
threshold parameter as the iteration goes on. Let us denote an image by /, its Fourier transform by F and partial 
Fourier samples by G. Here is the flow-chart of a generic image reconstruction algorithm via alternating projections. 



Algorithm 1. Image Reconstruction via Alternating Projections 

Input: observation data G and sampling pattern S; 
Output: reconstructed image / 

• Initialization: obtain f^^^ by ad-hoc back-projection method; 

• Main loop: for A; = 0, 1, kmax, 

- Projection onto prior constraint set: apply nonlinear filter to f^'^^ (at the customer's 
choice); 

- Projection onto observation constraint set: F^^\m,n) = G{m,n) for 
{{m,n)\S{m,n) = 1}; 



Assuming the function of nonlinear filtering is available as a module, the above algorithm can be easily 
implemented (fewer than ten lines of MATLAB codes). In our implementatior ^ we have reused the tools of PM 
diffusion, nonlinear diffusion, TI thresholding and SA-DCT thresholding. Fig. [T] includes the evolution of PSNR 
results for four different nonlinear filters. It can be observed that all of them can achieve perfect reconstruction 
(PSNR > A8dB) - the performance achieved by llmagic ( http ://ww w. acm.caltech. edu/I I magic/p with identical 
experimental setting. The behavior between TV-diffusion and TI thresholding is similar, which supports their 
theoretic connection established in |[T2l . What is most interesting observation seems to be the behavior of Perona- 
Malik diffusion - after a slow start, it rapidly catches up and the PSNR soai^s into the range of over 80dB (lossless 
up to some rounding errors). A similar "phase-transition" phenomenon has been observed for BM3D-based image 
reconstruction [31] (BM3D is a more recently-developed and computationally demanding tool of nonlinear filtering). 

IV. Discussions 

So what have we learned? I would argue we - signal processing researchers - can learn three important lessons 
from the above experiment: 

^The deeper reason is only partially known to be connected with overcoming the difficulty with getting trapped by local optimum of a 
nonconvex function. 

^The source codes can be accessed at |http://www.csee.wvu.edu/~xinl/code/TV_recon.rarj 
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Fig. 1. PSNR Profiles of Alg. 1 with different projection operators: a) PM diffusion [4J; b) nonlinear diffusion [5J; c) TI thresholding (6); 
d) SA-DCT thresholding Q. 



• All roads lead to Rome. Equipped with powerful analytical skills, mathematicians can prove the conditions for 
perfect reconstruction of smooth signals; equipped with powerful computational resources, engineers can achieve 
the same objective by recycling some of old ideas or tools. Mentally reproducible objects such as mathematical 
theorems and experimentally reproducible objects such as computational algorithms are simply a matter of different 
taste. Some analytically difficult hard-bones such as Perona-MaUk diffusion continues to prevail and offer nice 
surprises in numerical experiments. 

• Understand the hmitation of mathematical models. Statistician George Box once said, "All models are wrong; 
some are useful." The usefulness of any mathematical model relies on how well it matches real-world data. For 
computer-generated images such as phantom, we might find several seemingly-distant models (e.g., diffusion vs. 
thresholding) are intrinsically connected. It can be further argued that such kind of data represent pathological 
cases which one do not find in the real world and therefore could mislead our effort on mathematical modeling. 
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Nevertheless, the model underlying the CS theory is identical to the TV model despite the algorithmic differences 
between llmagic and TV-diffusion. 

• Use computational tools wisely. Like any other tool invented by humans, ii -optimization is appropriate for 
certain types of problems. The practice of treating every engineering problem hke a nail and attempting to solve it 
with a hammer (Zi -optimization) has shown to be less fruitful. Of course, we should not blame inventors of those 
tools by their creation; but realize that it is our own responsibility to choose the right tool to use and use it wisely. 
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